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Abstract 

In this study, we investigate the dispersive properties of smoothed particle 
magnetohydrodynamics (SPM) in a strongly magnetized medium by using 
linear analysis. In modern SPM, a correction term proportional to the di¬ 
vergence of the magnetic fields is subtracted from the equation of motion 
to avoid a numerical instability arising in a strongly magnetized medium. 
From the linear analysis, it is found that SPM with the correction term suf¬ 
fer from significant dispersive errors, especially for slow waves propagating 
along magnetic fields. The phase velocity for all wave numbers is signifi¬ 
cantly larger than the exact solution and has a peak at a finite wavenumber. 
These excessively large dispersive errors occur because magnetic helds con¬ 
tribute an unphysical repulsive force along magnetic helds. The dispersive 
errors cannot be reduced, even with a larger smoothing length and smoother 
kernel functions such as the Gaussian or quintic spline kernels. We perform 
the linear analysis for this problem and hnd that the dispersive errors can 
be removed completely while keeping SPM stable if the correction term is 
reduced by half. These hndings are conhrmed by several simple numerical 
experiments. 
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1. Introduction 


Smoothed particle hydrodynamics (SPH) is an entirely Lagrangian parti 
cle method for simulating fluid flows [l|,l2j. This Lagrangian nature offers 


ma¬ 


jor advantages when SPH is applied to problems with a large, dynamic range 
of spatial scales. Furthermore, SPH can easily incorporate other physics such 
as self-gravity, radiative transfer, or chemistry. Thus, SPH is widely used in a 
variety of astrophysical problems such as formation of large-scale structures, 
galaxies, stars, and planets. 

Recently, several authors have tried to extend SPH to magnetohydrody¬ 
namics (MHD) because magnetic helds play an important role in a variety of 
astrophysical environments. In this study, we call SPH for MHD “smoothed 
particle magnetohydrodynamics” (SPM). Price and Monaghan has de¬ 
veloped an SPM method with artihcial viscosity and resistivity [also see 1^ . 
Iwasaki and Inutsuka have applied Godunov’s method to SPM. We call 
it “GSPM”. Recently, Iwasaki and Inutsuka have modihed their original 
GSPM formulation, based on their derivation of the equation of motion in 
GSPM from an action principle. 

Unfortunately, conservative formulations of SPM are known to inevitably 
suffer from numerical instability for low (3 plasma because of negative stress, 
where [3 is the ratio of the gas pressure to the magnetic pressure. This insta¬ 
bility has already been pointed out by Phillips and Monaghan [also see 
[^. Among several methods proposed for removing the numerical instability 
n, S, [^ , an approach by Bpiw^ Omang, and Trulsen is widely used in 

loll. A broad discussion of stable SPM 


modern SPM methods [e.g., i, 
schemes is found in a review by Price 


In the approach of Bprve, Omang, 
and Trulsen [9t|, a correction term, B(V -Bj/dvr, is subtracted from the right- 
hand side of the equation of motion. The correction term is essentially zero 
if V ■ B = 0 is satished. By a linear analysis of SPM, Bprve, Omang, and 
Trulsen [l^ (hereafter BOT04) found that half of the correction term, or 
B(V • Bj/Svr, is large enough to remove the numerical instability. This was 


conhrmed by Barnes, Kawata, and Wu [13|, who found that half the correc¬ 


tion term provided the least error and minimized the violation of energy and 
momentum conservation in a variety of test calculations. Bprve, Omang, and 
Trulsen [l^ have proposed a sophisticated method, wherein the size of the 
correction term varies among the SPH particles. 

However, it is still unclear how the correction term affects the capability 
of SPM to accurately model fluid flows. Guiding the optimal selection of the 
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151 have 


amount of correction in a rigorous manner is important. Balsara 
investigated the linear stability of various SPH formulations, kernel functions, 
and ratios of smoothing length to interparticle distance [also see B0. They 


have suggested an optimal range of parameters from obtained dispersion 
relations. Although these analyses are valid only for the linear regime and 
a regular particle conhguration, they provide us with a lot of knowledge for 
achieving further improvements to SPM schemes. Pioneering work for SPM 
has been done by BOT04, who investigate the dispersive properties of SPM. 
They parameterize the size of the correction term with a free parameter ^ 
(0 < ^ < 1), and use ^ x BV ■ B/dvr as the correction term. They found 
that as mentioned above, ^ = 1/2 is large enough to stabilize SPM for low 
(3 plasma, and also suggested that smoother kernels, such as Gaussian or 
the quintic spline kernels, reproduce correct phase velocities, while the cubic 
spline causes large dispersion errors. However, their study is restricted to 
several linear waves in the long-wavelength limit and to the case with .^ = 1/2 
although many authors still adopt = 1- 

In this study, we investigate detailed dispersive properties of SPM for 
low /3 plasma by changing the kernel functions, and the ratio of smoothing 
length to interparticle distance. From the results, a suggestion of an optimal 
choice of the size of the correction term is provided. 

The paper is organized as follows: in Section |2l the basic equations of 
SPM are reviewed. In Section [31 a dispersion relation is derived from the 
basic equations of SPM and its asymptotic behavior in the long- and short- 
wavelength limits is discussed. The results of the linear analysis are presented 
in Section 01 To conhrm the results of the linear analysis, several numerical 
experiments are demonstrated in Section [5l Our results are discussed in 
Section [HI Section [7| offers a summary. 


2. SPM Equations 

The basic equations of MHD are given by 

| + V.(pv) = 0. 

dv^ 1 f 

= — 


dt p 




and 


^ (t) ^ 


( 1 ) 

( 2 ) 

( 3 ) 


3 







where is the stress tensor, 


( R2\ 

and d/dt = d/dt + v - V is the Lagrangian time derivative. The second term 
on the right-hand side of eqnation ([2]) is the correction term introdnced to 
remove the nnmerical instability (BOT04). The parameter ^ specihes the 
size of the correction term. In this stndy, ^ is assnmed to be constant for 
all particles. For simplicity, instead of the energy eqnation, the isothermal 
eqnation of state is assnmed: 


P 


c^p, 


( 5 ) 


where c is the sonnd speed. In the adiabatic case, the dispersive properties 
of SPM are expected to be qnalitatively the same as those in the isothermal 
case. 

In SPH, the density of the a-th particle is evalnated by the following 
eqnation: 

pa = '^mbW{xa-^b,h), ( 6 ) 

b 

where the snbscripts denote the particle label, mb is the mass of the &-th 
particle, W (x, h) is a kernel fnnction, and h is the smoothing length. In the 
linear analysis presented in this stndy, the smoothing length is assnmed to 
be constant. In the nnmerical experiments shown in Section 0 a variable 
smoothing length is nsed. 

There are several conservative formnlations of SPM. Here, we show two 
schemes: the standard SPM formnlation, the GSPM formnlation. The basic 
eqnations of standard SPM BSB are given by 


dt 


rjifiu \ 


and 

d_ 

dt \ p J 


± El] 

dt\p)a ^hpl 




( 8 ) 
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aJ diss 












where Wab{,ha) = hh(xa — :x.b,ha), the artificial dissipation terms (viscosity 
and resistivity) are denoted by the subscript of “diss”, and hi denotes the 
effect of the variation of the smoothing length. The detailed expression is 
described in Price and Monahan j^. 

Iwasaki and Inutsuka [^, implemented GSPM given by 

^ - Vm (SiLv-W (h 1 +fidv-M/ (h d 

- + ( 9 ) 


in which the quantities with the asterisks are the results of a Riemann solver. 
The evolution equation of /p is the same as that of standard SPM, except 
for the artificial resistivity term. In Iwasaki and Inutsuka j^, the resistivity 


term is derived from the method of characteristics for Alfven waves [17 
Note that equations ([9]) is a simpli^e^ version of a rigorous expression which 
is found in Iwasaki and Inutsuka 


Balsara 0 and Morris ai 0 have shown that the choice of kernel func¬ 
tions significantly influences on dispersive properties of SPH, and should be 
determined by the requirements of accuracy and computational efficiency. 
We consider several kernel functions to investigate how they affect the dis¬ 
persive properties of SPM. The first is the Gaussian kernel, which has the 
best interpolation accuracy. It is given by 


WG{r, h) = 


y/irh) 


\ 


(r/hp 


( 10 ) 


where d is the number of dimensions. One disadvantage of the Gaussian 
kernel is its infinite extent; in actual calculations, the contributions for r > 
3.1h are ignored. Also, in the linear analysis, we use a truncated Gaussian 
kernel at r = 3.1h that is denoted by WtG- Most works use the cubic spline 
kernel (l^ . 1 ^ . which is given by 


W,{r, h) 


Cd 

¥ 


1 - 1^2 + |s3 
(2 - .)^ /4 
0 


for 0 < s < 1 
for 1 < s < 2 
otherwise 


( 11 ) 


where s = rjh^ and Ci = 2/3, C 2 = IO/Ttt, and G 3 = l/vr. This kernel 
has a great computational advantage, as it has a compact support at r = 
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2h. However, the interpolation is inaccurate and leads to relatively large 
dispersion errors in sound waves 
kernel, 


15[ |. We also consider the quintic spline 


W,{r, h) 


Ca 

¥ 


- 10 s® + 30s^ - 60 s2 + 66 
5s® - 45s^ + 150 s3 - 120s 2 + 75s + 51 
-s® + 15s^ - 90 s3 + 270s 2 - 405s + 243 
0 


for 0 < s < 1 
for 1 < s < 2 
for 2 < s < 3 ’ 
otherwise 


( 12 ) 

where Ci = 1/120, C 2 = 7 / 4787 r, and ^3 = l/1207r. This kernel is smoother 
and more accurate kernel than IT 4 . 


3. Linear Analysis 

A two-dimensional (2D) rectangular lattice of particles with an interval of 
Ax is considered as an unperturbed state. The position of the a-th particle 
position is given by x^o = (Ax) a, where a is the 2D integer vector, a = 
(oa;, tty) {ttx, tty = 0,1, 2, • • •). The particle mass mo and smoothing length h 
are assumed to be constant. The following perturbations are considered: 


Xa = Xao + 5Xa 

(13) 

Pa Po T ^Pa 

(14) 

Va = 

(15) 

Ba = Bq -|- 5B(i 

(16) 

Pa = Po + C^Spa, 

(17) 


where the subscript “ 0 ” indicates physical variables in the unperturbed state. 
For simplicity, the unperturbed magnetic held is assumed to be parallel to 
the x-direction, and huctuations in the x-direction that correspond to Alfven 
wave are not considered, although their propagation can be partly inves¬ 
tigated by fast waves propagating along magnetic helds for low (3 plasma. 
This means that we consider fast and slow modes that oscillate in the (x, y)- 
plane. We assume that the perturbations have the following space and time 
dependence: 

OC exp {i (k ■ Xqo - Wt)} , (18) 

where 6Qa = (^x^, (fv^, ^B^). 
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Substituting equations (lT^ - flT7|) into the standard SPM equations (E]), 
m, (0), and dyia/dt = ^a, we obtain a dispersion relation. The artihcial 
dissipation term is omitted in this analysis. We do not repeat the detailed 
derivation of the dispersion relation that was already shown by BOT04. The 
dispersion relation is given by 

det + A^'') = 0 (19) 


where “det” indicates determinant, 

Po Po I 47r J 

( 20 ) 

is the unperturbed stress tensor modihed by the correction term, 


rjnfJ.h' _ 


- P< 




( 21 ) 


0 ^^, and are given by 


^ Po 


— e“**^4XaO-Xi,o)j 


dW, 


ab,0 


dKo 


( 22 ) 


and 


ht^ — 


Ef (1 + 

K Po 


,-ik-(Xa0-X60 


') 


SKo 




^-ik-(xao-X6o 


Po 


>) 


SKijSxm' 


(23) 


(24) 


The hrst term on the right-hand side of equation fl20l) comes from the per¬ 
turbation of (1/Pa + 1 /Pfe)and the second term comes from the per¬ 
turbation of in equation (17|) . 

Note that the dispersion relation derived from the linearized GSPM equa¬ 
tions is identical to equation flT^ if the physical quantities with an asterisk 
are evaluated at the arithmetic mean between the a- and 6-th particles. Thus, 
this linear analysis is valid both in standard SPM and GSPM. 

Since equation flT^ is a quadratic equation in terms of two modes are 
obtained. In this study, the mode with larger (smaller) cj^ is referred to as 
the fast (slow) mode. The numerical phase velocities Ci;/|k| of the fast and 
slow modes are denoted by Cf,num(k) and Cs,num(k), respectively. 
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3.1. Asymptotic behaviors 

In this section, we investigate the asymptotic behavior of the dispersion 
relation in the short- and long-wavelength limits. 


3.1.1. Short-wavelength Limit 

SPM without the correction term is unstable for low /plasma because 
negative diagonal components appear in the stress tensor [sj]. The develop¬ 
ment of numerical instability begins with the growth of fluctuations whose 
scales are comparable to Ax. Thus, in this section, we investigate the asymp¬ 
totic behavior in the short-wavelength limit and the size of ^ required to 
stabilize SPM. Note that this already has been done by BOT04, who numer¬ 
ically solve the dispersion relation. Here, we show that their conclusion is 
reproduced in the following simple analytical manner. 

In the discrete system, the largest wavenumber is fc = tt/Ax where the 
unit wavelength is expressed by two particles. A compressible wave propa¬ 
gating along Bo (k = (tt/Ax, 0)) is considered. This corresponds to the slow 
mode. In this case, we obtain 0^ = -0^ = 0 from equations fl22|) and ([23]). 
From equation (fUj), the dispersion relation becomes 



where (3 = SttPq/Bq. The summation in equation fl25p is positive in normal 
situations j8|. Without the correction term = 0), one can see that SPM 
becomes unstable, since < 0 if /5 < 1. To be stable, the following condition 
should be satished: 

e>Un=^(l-/9), (26) 

where .^min is the minimum value of ^ needed to ensure stability for a given 
/3. This linear dependence of ^min on f3 is the same as that found numerically 
in BOT04 (see their Fig. 7). In the strong magnetic held limit (/3 —)• 0), 
^ = 1/2 is large enough to stabilize SPM. 


3.1.2. Long-wavelength Limit 

In the long-wavelength limit (|k|Ax 0), summations can be replaced 
by integrals in equation 0191) 20| . For instance, is approximated by 






dW (x, h) 2 


dx^ 


d^x = -Tfc^IT(k), 


(27) 







where integration by parts and hh ^ 0 for a: —)■ oo are used, and f^(k) is 
the Fourier transform of hF(x, h), given by 



FF(x, h)(fx. 


(28) 


In the similar way, one obtains 


(j)^^ ~ ik^^W and ~ 


(29) 


Using equations (127|) and (12^ . equation (120|) becomes 



If the Gaussian kernel is applied, its Fourier transform is also Gaussian, 


IU(k, h) = e Thus, IF —)■ 1 + 0{{hk)‘^) for |k|h —)■ 0. Using this fact, 

the dispersion relation (equation (IT^ l becomes 



(31) 


This dispersion relation holds for both fast and slow waves. Note that equa¬ 
tion fl3Tll does not contain Thus, waves in the long-wavelength limit are 
not affected by the correction term. Ideally, the dispersion relation of SPM 
satishes the correct phase velocity as long as equations fl2711 and fl29ll are 
valid. 

Gomparing equations fl5U]) and (1201) . one can see that the correct phase 
velocities come from the second terms on the right-hand sides. The first 
term on the right-hand side of equation (l30|) becomes zero because IF ~ 1. 
In reality, the numerical dispersion relation (IT^ is expected to deviate from 
that in equation fl3ip because finite discretization errors are introduced in 
equations flTT)) and fl2^ . The first term on the right-hand side of equation 
fl20|) causes larger dispersive errors than the second term. This is because 
contains the second derivative of the kernel function that has larger errors 
(see equation flMjl b In particular, the cubic spline kernel leads to large errors 
in because its second derivative is a broken line. Thus, dispersive errors 
mainly come from the hrst term on the right-hand side of equation fl2U]) . 
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4. Results 


4 . 1 . Phase Velocities in Long-wavelength Limit 

In this section, the dispersion relation fll9p is solved by considering a 
sufficiently small wavenumber (|k| = 10 “^ 7 r/Ax) and changing the angle 6 
between k and Bq in order to investigate whether SPM correctly reproduces 
the phase velocities shown in Section [3.1.21 

4 . 1 . 1 . Fast Mode 

Fig. [T] shows the numerical phase velocities of the fast mode as a function 
of the angle 6 for various values of /3, h, and various kernel functions. The 
value of ^ is assumed to be 1. The exact solutions are plotted by the gray 
lines. For low /3 plasma, the fast wave at 6 * = 0 represents the incompressible 
pseudo-Alfven wave whose phase velocity is Ca = Bq/ yjAn pq. The phase 
velocity increases with 9 and gradually changes into a compressible mode. 
At 9 = 7 r/ 2 , the phase velocity reaches a maximum value of + c^. In 
the upper column of Fig. [ 1 ], for h = Ax, all kernel functions reproduce the 
exact solutions within an error of 1% although the Gaussian kernel provides 
slightly worse results. For h = 1.2Ax, the results for the Gaussian kernel 
are almost identical to the exact solution (see the lower column of Fig. [T]). 
Only the results obtained with IF 4 suffer from relatively large errors. The 
results for ^ = 1/2 are not shown, because it is con&med that the numerical 
dispersion relations do not much depend on 

It is found that SPM can reproduce the phase velocity of the fast mode 
reasonably well regardless of This behavior can be qualitatively understood 
from the one-dimensional dispersion relation for kx = ^ {9 = n/2) given by 

= (c^ -I- 2cl) (T™ — (jp'ipy) -t- (c^ -f cl) . (32) 


As mentioned in Section [3.1.21 the correct phase velocity comes 

from the second term on the right-hand side of equation fl32|) and the dis¬ 
persion errors arise mainly from the hrst term. To obtain the correct phase 
velocity, the hrst term on the right-hand side of equation (13^ should be 
negligible compared with the second term. We can see that the coefficient 
of (^ra — (^Vipy) in equation fl3^ is comparable to that of cjp'ipy. Thus, as 
long as \ /{(IP'ipy)\ -C l is satished, the numerical phase velocities 

agree with the exact values within sufficiently small errors. For 9 7 ^ 7 r/ 2 , a 
term proportional to ^ appears in the hrst term on the right-hand side of 


10 









equation (15^ . Also in this case, the effect of the first term can be neglected, 
as the phase velocity of the fast mode is comparable to Ca- That is why the 
numerical phase velocities agree with the exact phase velocities reasonably 
well and do not depend much on ^ for all 9. 

4 . 1 . 2 . Slow Mode 

The slow mode exhibits completely different dispersion properties than 
the fast mode. Fig. [2] shows the numerical phase velocities of the slow mode 
as a function of the angle 6 for various values of fd and h, and various kernel 
functions, with ^ = 1. At 6* = 0, the slow mode corresponds to the sound wave 
propagating in the direction parallel to Bq. The phase velocity decreases with 
9 and gradually becomes the incompressible mode. At 6 * = 7 r/ 2 , the phase 
velocity becomes zero. First, we focus on the cases with h = Ax shown in the 
upper panels of Fig. [2l The results with IT 4 are significantly underestimated 
for fd = 0.1 and leads to a numerical instability ai jd = 0.01. Even for the 
smoother kernels (PFg, khtc, fFe), the dispersion errors are significantly large 
and becomes worse for smaller fd. Although only Wq provides a better result 
for fd = 0 . 1 , the stronger magnetic field [fd = 0 . 01 ) makes the results worse. 
From the lower left panel in Fig. |2l for fd = 0.1, larger h {h = 1.2Aa:) 
improves the results with the smoother kernels while the result with W 4 does 
not improve. Also, in the case with larger h, for the stronger magnetic field 
{fd = 0.01), the dispersion errors become worse except for Wq. In summary, 
the dispersion errors of the slow mode are large and becomes significant as fd 
decreases. A larger smoothing length makes the errors lower with every kernel 
except for the cubic spline, although large errors still remain for sufficiently 
low fd. Thus, larger smoothing length cannot be an ultimate solution. 

Fig. [3] is the same as Fig. [2] but setting ^ = 1/2. All kernels can 
provide the correct phase velocities for all cases although the results with IF 4 
have relatively large errors. BOT04 have already investigated the dispersion 
properties for ^ = 1 / 2 , and their results are consistent with ours. 

Figs. |2] and |3] reveal that SPM with ^ = 1 exhibits significantly large 
dispersive errors that disappear if ^ = 1/2 is used. This behavior can be 
qualitatively understood from the simple one-dimensional dispersion relation 
for ky = 0 {9 = 0) given by 

a;2 = {2c^ + cl{2( - 1)} (T"" - (33) 

To obtain the correct phase velocity, c, the first term on the right-hand side 
of equation fl32]l should be negligible when compared with the second term. 
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Figure 1: Phase velocity of the fast mode versus the angle between k and Bq for (/3, h) = 
(0.1, Ax), (0.01, Ax), (0.1, 1.2Ax), and (0.01, 1.2Ax). The results with ^ = 1 and ^ = 1/2 
are almost identical, and so only those with ^ = 1 are presented. The red, green, blue, 
and orange lines correspond to the results using the Gaussian, truncated Gaussian, cubic 
spline, and quintic spline kernels. In each panel, the exact solutions are shown by the gray 
line. 
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Figure 2: Phase velocity of the slow mode versus the angle between k and Bq for 
{l3,h) = (0.1, Ax), (0.01, Ax), (0.1,1.2Ax), and (0.01,1.2Ax). The red, green, blue, and 
orange lines correspond to the results using the Gaussian, truncated Gaussian, cubic spline, 
and quintic spline kernels. In each panel, the exact solutions are shown by the gray line. 
All results are taken at ^ = 1. 
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Figure 3: Phase velocity of the slow mode versus the angle between k and Bq for 
(/3, h) = (0.1, Ax), (0.01, Ax), (0.1,1.2Ax), and (0.01,1.2Ax). The red, green, blue, and 
orange lines correspond to the results using the Gaussian, truncated Gaussian, cubic spline, 
and quintic spline kernels. In each panel, the exact solutions are shown by the gray line. 
All results are taken at ^ = 1/2. 
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(a) ^=1.0 (fast mode) (b) ^=0.5 (fast mode) 



Figure 4: Color maps of the phase velocities of the fast mode in the {kx^x/TT, kyAx/ir) 
plane for (a)^ = 1 and (b)^ = 1/2. The Gaussian kernel is considered. The parameters /3 = 
0.1 and h = 1.2Aa: are used. In each panel, the color shows the numerical phase velocity 
normalized by the exact phase velocity depending on 9. The gray region corresponds to 
the unstable region. 

However, in contrast to the fast mode case, the coefficient of 
in equation is much larger than that of for low (3 (ca c). Thus, 
the hrst term can be important even if — (/)®'^^)/(0®'^*)| is sufficiently 
small. Note that the contribution of the magnetic held disappears only at 
^ = 1/2 in the hrst term of equation (jSS]), and the dispersion relation is 
reduced to that without a magnetic held. That is why the errors are mostly 
eliminated for ^ = 1/2, as shown in Fig. |3l 

These hndings in the long-wavelength limit suggest that the best choice 
for is 1/2, as SPM with .^ = 1 suhers from signihcant errors. 

4-2. Dispersion Relation 

In this section, the overall dispersive properties of SPM are investigated. 
Figs. and |Dd show color maps of the numerical phase velocities of the 
fast mode in the {kxAx/7T, kyAx/Ti) plane for ^ = 1 and 1/2, respectively. 
The smoothing length is assumed to be h = 1.2Ax, and the plasma /3 value 
is hxed to be // = 0.1. Here, we focus on the results with the Gaussian 
kernel. It is conhrmed that the qualitative properties do not depend on the 
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(a) ^=1.0 (slow mode) 
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(b) ^=0.5 (slow mode) 



0.0 0.5 1.0 
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Figure 5: Color maps of the phase velocities of the slow mode in the {kxAx/ir, kyAxIn) 
plane for (a)^ = 1 and (b)^ = 1/2. The Gaussian kernel is considered. The parameters P = 
0.1 and h = 1.2Aa; are used. In each panel, the color shows the numerical phase velocity 
normalized by the exact phase velocity depending on 9. The gray region corresponds to 
the unstable region. 

kernel functions. In Fig. IH the numerical phase velocities are divided by 
the corresponding exact solutions depending on 6 . In both panels in Fig. 01 
Cf,num/cf(^) has peaks at the origin and monotonically decreases with |k| at 
all 6 , where Cf is the exact phase velocity of the fast mode. The difference 
between the results with = 1 and 1/2 is found only in the region where 
|k|Aa;/7r > 0.5 and 6^ ~ 0. For p = 1/2, Cf,num/cf declines more rapidly than 
it does for = 1. This behavior will be explained later. Thus, the dispersion 
relation of the fast mode does not depend much on the value of p. This is 
consistent with the hndings in Section 14.1.11 

Next, the dispersion relations of the slow mode are investigated. Fig. |5] 
is the same as Fig. 01 but for the slow mode. First, the result with ^ = 1, 
as shown in Fig. |5k, is investigated. Fig. |5k shows that Cs^num/c exhibits 
anomalous features, although Cs,num can reproduce the correct phase velocities 
in the long-wavelength limits for the Gaussian kernel (see the lower left panel 
in Fig. 12]). Cs,num/cs increases from the origin with |k| especially in the 
direction parallel to Bq {0 ~ 0). Around |k|Aa;/7r ~ 0.4, Cs,num/cs reaches 
a maximum whose value is ~ 2.5 times larger than the exact solution. For 
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larger |k| (|k|Ax/7r > 0.4), Cs^num decreases with |k|. 

These anomalous features of the dispersion relation for ^ = 1 completely 
disappear at = 1/2, as shown in Fig. |5 ]d where Cs,num/cs monotonically 
decreases with |k| from the center. 

Fig. |5] reveals that the phase velocities around 6 = 0 strongly depend 
on To see this more clearly, the phase velocities at 6 = 0 and 6 = 0.05 
are plotted in Fig. E] as a function of |k|. Fig. shows the results with 
^ = 1. First, we focus on the case where 6 = 0, shown by the solid lines. 
As mentioned above, the numerical phase velocities agree with the exact 
values for both the fast and slow modes. From k = 0, Cf^num decreases with 
k while Cs^num increases. At kAx/n ~ 0.4, the two branches have the same 
phase velocity. Beyond this point, it can be clearly seen that the fast mode 
branch smoothly connects with the slow mode branch. This happens simply 
because larger (smaller) is referred to as the fast (slow) mode in this 
study, regardless of its eigenfunction. The “real” fast (slow) mode should 
be an incompressible (compressible) mode at 6 = 0. By investigating the 
corresponding eigenfunctions, it is conhrmed that the “real” fast and slow 
modes intersect at kAx/n ~ 0.4. This means that the slow mode changes 
into an incompressible mode and the fast mode changes into a compressible 
mode beyond the intersection point. The dashed line in Fig. [6] indicates the 
results with 6 = 0.05. We can see the mode exchange between the fast and 
slow branches around kAx/ir ~ 0.4. From the eigenfunctions, it is conhrmed 
that the fast (slow) branch is changed into a compressible (incompressible) 
mode by the mode exchange. Fig. [6^ indicates that the phase velocity of 
the compressible mode (the “real” slow mode) is supersonic (> Cg) for all 
wavenumbers. 

Fig. [7] shows the /9-dependence of the maximum phase velocity of the 
compressible mode for .^ = 1 and 6 = 0. As shown in Fig. |H^, the “real” slow 
mode has an off-center peak. From Fig. [TJ we can see that the maximum 
phase velocity increases with decreasing f3. The maximum phase velocity for 
low {3 is well htted by Ca/2 shown by the dashed line Fig. [71 This dependence 
clearly indicates that the enhancement of the phase speed of the compressible 
wave comes from the magnetic helds. 

Fig. [6 ]d shows the results for ^ = 1/2. The phase velocity of the slow 
mode monotonically decreases with k. Also for .^ = 1/2, the mode exchange 
occurs at larger k. From Figs|6K andjUb, the extended feature around 6^0 
in Cf^num in Fig- SK corresponds to a compressible mode with a large phase 
velocity (> Cg) created by mode exchange. 
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Figure 6: Dispersion relations for (a)^ = 1 and (b)^ = 1/2. The Gaussian kernel is used. 
In each panel, the blue and red lines corresponds to the fast and slow modes, respectively. 
The solid and dashed lines indicate the results with ^ = 0 and 6 = 0.05, respectively. 
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Figure 7: The /3-dependence of the maximum phase velocity of the compressible mode 
for the ^ = 1 and 9 — 0 case. The Gaussian kernel is used. The dashed line indicates Ca/2. 
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Figure 8: Dispersion relations for (a)^ = 1 and (b)^ = 1/2. In each panel, the blue and 
red lines correspond to the fast and slow modes, respectively. The solid and dashed lines 
indicate the results with 6 = 0 and 6 = 0.05, respectively. 

5. Numerical Experiments 

Fig. [6] reveals that the dispersion relation is abnormal if ^ = 1 is used. 
In Section 14.21 it is found that the compressible waves at short-wavelength 
propagate with supersonic velocities. In this section, to test this dispersive 
property, three simple test calculations are demonstrated. 

5.1. Propagation of an Isolated Wave 

Fist, test of the propagation of an isolated wave is performed. This is 
a severe test for the propagation of linear waves because an isolated wave 
is composed of many linear waves with different wavelengths. Thus, if the 
sound speed numerically has a A;-dependence, the shape of an isolated wave is 
expected to change during its propagation. In other words, the deformation 
of an isolated wave clearly shows the dispersion errors. 

In the unperturbed state, the particles are distributed in a rectangular 
lattice in the domain of [—2,2] x [—1,1]. The density and pressure are as¬ 
sumed to be 1. A uniform magnetic held is introduced in the x-direction. 
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and its strength is set snch that (3 is 0.1. A velocity pertnrbation given by 
5v = 0.01 exp(—(a:/3h)^) is added. The two isolated waves propagate ont- 
ward from the origin. As the sonnd wave does not have dispersion, the two 
isolated waves shonld keep their shape with the sonnd speed (c = 1 ) as they 
propagate. Cha and Whitworth [2l| did the same test calcnlation withont 
magnetic helds. 

Fig. [H^ shows a snapshot of the velocity pertnrbation at f = 0.6 for = 1. 
Only the isolated wave propagating rightward is plotted. The gray lines 
indicate the exact solntion. The colors indicate the difference of the kernel 
functions. For ^ = 1, with all kernel fnnctions, the isolated waves break 
into smaller waves that propagates at snpersonic velocities. This dispersive 
property is consistent with the resnlts of linear analysis. The resnlts with 
Wq are almost identical to that with Wtc- The isolated wave with IT 4 shows 
a larger speed, as expected in Fig. [5l On the other hand, for = 1/2, there 
is no destrnction of the isolated wave, and the results agree with the exact 
solution quite well, although the result with shows a slightly larger phase 
velocity. This is explained by the linear analysis (see the lower left panel of 
Fig. [3]). These hndings suggest that the optimal choice for is 1/2. Otherwise 
SPM provides completely incorrect results for the wave propagation. 


5.2. Colliding Flow Test 

The linear analysis is valid only for the linear waves. In this test, we con¬ 
sider a colliding flow test that involves shock waves. We investigate whether 
the dispersive errors affect the shock structures or not. The computational 
domain is —1 < a; < 1, —0.0625 < y < 0.0625, and the periodic bound¬ 
ary condition is imposed in the ^/-direction. Initially, the density is uni¬ 
form, and the gas moves toward a; = 0 from both directions with a velocity 
±3.75. The corresponding Mach number is 4. The initial magnetic held is 

B = ^a/Stt//?, oj. As the magnetic held is parallel to the x-direction, it 
should not ahect the gas dynamics. The truncated Gaussian kernel is used. 
We consider two cases: (3 = 0.1 and 0.01. In the exact solution, two shock 
waves propagates outward from x = 0. In this test, we investigate how the 
strong parallel magnetic held numerically ahects the gas motion. The initial 
particle distribution is set to be a random distribution that is relaxed until 
the density dispersion is sufficiently small. In this test, a Riemann solver is 
used to capture shock waves. The hyperbolic divergence cleaning method 
(^. is also used. 
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Fig. inti- shows the results with ^ = 1 and ^=l/2for/3 = 0.1. The black 
line indicates the exact solution. Both cases agree with the exact solution 
quite well. This is because the shock jump condition is determined only 
by mass and momentum conservation laws, regardless of their dispersive 
properties. For the stronger magnetic held case {(3 = 0.01), the effect of 
the dispersive errors is shown in Fig. [Hb. For ^ = 1/2, the density of the 
shocked gas can reproduce the exact solution reasonably well within an error 
of 2%. This small error occurs due to the perpendicular magnetic held By 
numerically generated at the shock front. It works as an additional pressure, 
leading to the smaller density jump. For ^ = 1, the density prohle is quite 
diherent from the exact solution. The shock fronts broaden and small waves 
propagate toward the upstream with a supersonic velocity larger than the 
converging velocity 3.75. This can be understand by Fig. [71 For [3 = 0.01, 
the maximum phase velocity of the compressible wave is 7, which 

is larger than the converging velocity. Thus, waves can propagate toward 
upstream against the converging how. 

From this converging how test, it is found that given /?, there is a mini¬ 
mum shock speed below which waves with short wavelengths propagate up¬ 
stream and disturb the preshock gas. The minimum shock speed corresponds 
to the maximum phase velocity of the compressible wave, or Ca/2 derived from 
the linear analysis (see Fig. [7]). In other words, if the Alfven Mach number 
is smaller than 0.5, the dispersive errors are serious. 

5.3. Hydrostatic Equilibrium Under An External Gravity 

The dispersion errors originate from the fact that a parallel magnetic held 
numerically works as an additional repulsive force. This indicates that the er¬ 
rors can be important not only for waves but also for hydrostatic structures. 
Thus, in this test, we investigate whether SPM reproduces a hydrostatic 
structure under an external gravity, Qx = —2tanh(a;). As the initial condi¬ 
tion, we consider a uniform static gas with a uniform magnetic held in the 
a:-direction. The amplitude of the magnetic held is -^/Stt//?. The calculation 
domain is — 2 < a: < 2 and —1/4 < y < 1/4. The truncated Gaussian kernel 
is used. A periodic boundary condition is imposed in the ^/-direction, and the 
wall boundary condition is imposed in the x-direction. The initial particle 
conhguration is set to be a settled random distribution. Because of the exter¬ 
nal gravity, the plasma accumulates toward the midplane (x = 0). Finally, 
the density prohle is expected to relax to the hydrostatic equilibrium prohle, 
Peq{x,y) = 1/ cosh(a;)^. To avoid an undesirable oscillation in the relaxation 
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Figure 9: Density profiles at t = 0.2 for (a)/3 = 0.1 and (b)/3 = 0.01. Each point 
corresponds to an SPH particle. In each panel, the red and blue lines corresponds to the 
results with C = 1 and ^ = 1/2, respectively. The green lines indicate the exact solution. 


phase, we add a small drag force, —0.005v/At, to the equation of motion. 
If the maximum velocity becomes smaller than 0.001, the calculations are 
terminated. 

The upper panel of Fig. [TOh shows the obtained density prohles for 
(5 = 0.1. The lower panel of Fig. ITOk indicates the fractional residual, 
p{x) /peq^x) — 1. The result with ^ = 1/2 agrees with the hydrostatic prohle 
within sufficiently small error. The small deviation around x = ±2 comes 
from the boundary condition. Also for ^ = 1, SPM reproduces the hydro¬ 
static prohle reasonably well, although the density in the central region is 
underestimated and the low density tails are overestimated. However, this 
tendency becomes prominent for (3 = 0.01, as shown in Fig. ITOb . The density 
prohle with = 1 exhibits a more extended prohle than peq(|/)- On the other 
hand, even for /3 = 0.01, SPM with .^ = 1/2 can produce the correct prohle. 
From the test, it is found that the hydrostatic prohles along the magnetic 
helds are broadened by the artihcial repulsive force if .^ = 1 is used. 
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Figure 10: (Upper panels) Density profiles for (a)/3 = 0.1 and (b)/l = 0.01. (Lower 
panels) Fractional residual, p(x)/peq — 1. Each point corresponds to an SPH particle. In 
each panel, the red and blue lines corresponds to the results with ^ = 1 and ( = 1/2, 
respectively. The green lines indicate the hydrostatic profile. 
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6. Discussion 


6.1. Stability Against Particle Disorder 

The test calculations demonstrated in Section O show that SPM with 
^ = 1/2 removes dispersion errors. However, in a realistic situation where a 


blast wave propagates in a strongly magnetized medium, Tricco and Price [10 


(hereafter TP12) found that SPM with .^ < 1 produces numerical fluctuations 
behind the slow shocks and at the contact discontinuity. They found that a 
value of .^ = 1 leads to stable results. This is because the repulsive force is 
weaker for SPM with smaller ^ , leading to the particle disorder. 

To investigate the stability against the particle disorder in SPM with 
^ = 1/2, we perform the same blast wave test as in TP12. The total particle 
number is 512 x 512. Fig. ITTh and ITTb show the density maps at t = 0.03 
for ^ = 1 and ^ = 1/2, respectively. The truncated Gaussian kernel is used. 
One can see that the result with ^ = 1/2 is quite similar to that with ^ = 1. 
As shown in Section 15.21 SPM can treat shock waves even for ^ = 1 as long 
as the Alfven Mach number is larger than 0.5. Thus, it is difficult to identify 
the dispersion errors in this test. 

Note that the signihcant particle disorder found by TP12 does not appear 
in Fig. ITTb . Although this discrepancy may come from the difference of 
treatment of the numerical dissipations between GSPM and their scheme, 
the exact reason is still uncertain. However, also in the GSPM, the result 
with ^ = 1 is smoother than that with ^ = 1/2. Thus, it is possible that 
GSPM with ^ = 1/2 suffers from the serious particle disorder in more extreme 
situations. The best choice of ^ depends on situations from the point of view 
of accuracy and stability. As shown in the linear analysis, the dispersive 
errors are serious in slow waves propagating along magnetic helds. Thus, 
for example, to simulate a sub-Alfvenic turbulence in low fd plasma, a value 
of = 1/2 should be adopted. On the other hand, in dynamical situations 
where strong shock waves are important such as the blast wave test, the 
dispersive errors are not serious if the Alfven Mach number is larger than 0.5 
(see Section lA^ . In these cases, a value of ^ = 1 is acceptable. 


6.2. Comparison with Other SPM Formulation 

Besides the approach by Bprve, Omang, and Trulsen j^, Morris 16 
posed the following formulation; 


pro- 


dt 


-E 


mb 


P _L 


^aPl 


’-VaWabiha) + 


p j_P p2 

' '’■V'^Wabihb) 


^bPl 


24 












Figure 11: Density maps of the blast wave test at t = 0.03 for (a)^ = 1 and (b)^ = 1/2. 
The truncated Gaussian kernel is used. A value of ^ = 1/2 is used. 

+ 1 W ( BjBi - + v-Mh) 

4ir “ * V Po.Pi } 2 

In his formulation, the conservative form is used in the isotropic part of 
the stress tensor while the non-conservative form is used in the remaining 
part. Thanks to the non-conservative term, his formulation is free from the 
numerical instability. In this section, we compare the dispersion relations 
between the Morris formulation (SPMMoms) and the conservative form with 
the correction term (SPMcorr)- Linearizing equation fl34p . one obtains the 
following dispersion relation; 

det = 0, (35) 


where 




^ - Bsr) . 


0 ^ 

Po 


For + T [lifpr 
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(36) 


25 


density 



















Morris formulation 



k^Ax/jc 


Figure 12: Color maps of the phase velocities of the slow mode using SPMMorris in the 
{kx^x/TT,kyAx/TT) plane. The Gaussian kernel is considered. The parameters /3 = 0.1 
and h = 1.2Aa: are used. The color shows the numerical phase velocity normalized by the 
exact phase velocity depending on 9. 

Fig. [12] shows the color map of the numerical phase velocities of the slow 
mode using SPMMorris in the {kxAx/n, kyAx/n) plane. The smoothing length 
is assumed to be h = 1.2Ax, and the plasma jS value is hxed to be /3 = 
0.1. The Gaussian kernel is considered. One can see that the behavior of 
Cs,num is Quite similar to that in SPMcorr with = 1 (see Fig. [5k). The 
dispersion relations have an off-center peak of Cg^num around 0 ~ 0. Thus, the 
SPM formulation by Morris jl^ also suffers from the dispersive errors. This 
behavior can be qualitatively understood from equation (136|1 . The hrst term 
on the right-hand side of equation fl36|) that leads to most of the dispersive 
errors becomes large compared with other terms for low /3 plasma. Note that 
in SPMcorr fhe corresponding term proportional to is negligible 

if .^ = 1/2 is used. 

7. Summary 

In this study, we have investigated the dispersive properties of SPM with 
a correction term introduced to remove numerical instability in a strongly 
magnetized medium j^. The size of the correction term is parametrized by 
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^ (see equation ([2])). The findings in this study are summarized as follows: 

• As numerically found by Bqrve, Omang, and Trulsen jl^, the minimum 
value of ^ for removing the numerical instability is analytically derived 
as a function of the plasma (3 value in Section I3.1.1I 

• For the fast modes, it is found that SPM can reproduce correct phase 
velocities regardless of The dispersion properties are similar to those 
without magnetic fields. 

• The phase velocities of the slow modes is shown to significantly de¬ 

pend on For ^ = 1, which is used in most schemes, it is found that 
SPM suffers from significant dispersion errors with all kernel functions. 
The dispersion errors become worse for a lower value of f3. In the 
long-wavelength limit = lO^Ax/vr), the numerical phase veloci¬ 

ties are largely different from the exact values especially for the cubic 
spline kernel. A larger smoothing length and a smoother kernel func¬ 
tions are not ultimate solution if (3 is sufficiently small (see Fig. ED. 
The dispersion relations have an off-center peak of the phase velocities 
around \k\Ax/TT ~ 0.4 and 6* ~ 0, where 9 is the angle between k and 
Bq. Furthermore, the phase velocities are supersonic for all wavenum¬ 
bers above the Nyquist wavenumber tt/Ax (see Fig. [5D. The reason for 
this anomalous behavior is that a parallel magnetic field numerically 
works as an additional repulsive force and the phase velocities of the 
slow wave become supersonic. This fact can be understood by examin¬ 
ing the dispersion relation with respect to the slow wave propagating 
along a magnetic field. 

On the other hand, for ^ = 1/2, the dispersion errors found for ^ = 1 
completely disappear. This can be understood analytically from the 
dispersion relation. 

• To confirm the findings of the linear analysis, several simple numerical 
experiments are demonstrated. For the tests of linear isolated waves 
and hydrostatic equilibrium, SPM with ^ = 1 leads to unphysical re¬ 
sults while the exact solutions can be reproduced for .^ = 1/2. On the 
other hand, SPM with = 1 can treat the parallel shock waves if Alfven 
Mach number is larger than 0.5. This is because the shock condition is 
determined only by the conservation properties and slow waves cannot 
propagate toward upstream. If the shock speeds are smaller, waves 
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propagate upstream against the flow and cause disturbances. These 
results are consistent with the linear analysis. 

This study clearly shows that corrected SPM with ^ = 1 over-stabilizes 
the numerical instability and signiflcantly modifles the dispersion properties 
of the slow modes. To eliminate such dispersion errors, we suggests that 
^ = 1/2 is the best choice. These abnormal dispersion properties have not 
been found in the previous works, (e.g., the blast wave tests used to test the 
capability of schemes for low (3) because, as shown in Section 15.21 SPM can 
treat shock waves even for .^ = 1 as long as the shock speed is large. Thus, it 
is difficult to identify the dispersion errors found in this study from the blast 
wave test. The errors can be serious; for instance, in sub-Alfvenic turbulence 
for low (3 plasma, the phase velocity of the waves propagating along magnetic 
fields may be significantly overestimated. 

The linear analysis and the numerical experiments suggest that the best 
choice of ^ is 1/2 from the point of view of accuracy. However, as discussed in 
Section l6Tl SPM with ^ = 1/2 tends to less stable against particle disorder 
because of the small repulsive force. In dynamical environments where strong 
shock waves are important, a value of ^ = 1 is acceptable if the Alfven Mach 
number is larger than 0.5. 
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